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A Locally Implicit Method 
For Fluid Flow Problems 

1. Introduction 


The fluid flow inside the space shuttle main engine (SSME) traverses through a com- 
plex geometrical configuration. The flow is compressible, viscous and turbulent with pock- 
ets of separated regions. Computation of such a flow field requires the application of the 
state-of-the-art CFD technology. Several computer codes^ 1-4 ^ are being developed to solve 
three dimensional Navier-Stokes equations with different turbulence models for analyzing 
the SSME internal flow. Most of these computer codes require extensive computational 
effort and are rather slow to converge if accurate solutions are to be obtained in three 
dimensional domains. This report outlines an algorithm which has the potential of solving 
fluid flow problems with rapid convergence to steady-state solutions. 

The computational methods used in the Navier-Stokes codes fall into two major cat- 
egories: finite difference and finite element methods. Each have their own strengths and 
weaknesses. Some of the algorithms are designed to solve the unsteady Compressible 
Navier-Stokes equations, either by explicit or by implicit factorization methods, for sev- 
eral hundred or thousands of time steps to reach a steady-state solution asymptotically. 
Other algorithms attempt to reach the steady-state solution by relaxation methods. All of 
them require body-fitting curvilinear grids with sufficient resolution. Grid requirements, 
however, differ greatly with the algorithm used. Implicit factorization methods, typically 
use numerical transformations whereby the transformed grid in the computational space is 
uniform and rectilinear. This requires the grid to have indices which are separable in three 
directions for three dimensional problems, and also be reasonably smooth. However, such 
requirements may introduce grid singularities when complicated domains are discretized. 
Flow solver algorithm will have to deal with such grid singularities. Explicit schemes and 
finite element algorithms have less stringent requirements on the grid structure. However, 
their efficiency and accuracy also depend heavily on the grid structure. 

Explicit algorithms for solving time-dependent Navier-Stokes equations are usually 
bound by stability restrictions and take thousands of time steps to converge to a steady- 
state for three-dimensional flow' problems. Most of the current implicit algorithms use 
approximate factorization techniques and require index separable grid structure. The 
approximate factorization techniques are reasonably efficient for two-dimensional problems, 
but can be quite slow in convergence for three dimensional problems. We propose a locally 
implicit algorithm, which is based on a relaxation procedure, for solving time dependent 
partial differential equations to obtain steady-state solutions. The basic algorithm which 
is locally implicit but globally explicit is computationally efficient and is outlined in the 
next section. It is suitable for the application of convergence acceleration procedures such 
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as the multi-grid method which is also outlined. 

Test calculations with the current method for solving an elliptic partial differential 
equation, the diffusion equation, which simulates many features of the subsonic flow prob- 
lems are presented in the third section. Both single grid and multi-grid results are pre- 
sented. Solutions for inviscid and viscous flow problems from subsonic to supersonic flow 
problems are being computed successfully by this method and those results will be pub- 
lished in due course. 


2. Locally Implicit Scheme 

A relaxation scheme known as the locally implicit scheme is described below. The 
compressible Navier-Stokes equations can be written in the conservation form. 

dG 

aT +v/=0 < 21 > 

where G is the conservation variable vector and F = ( F x , F y , F z ) represents the flux vectors, 
including inviscid and viscous terms. Addition of the equation of state makes this system 
complete, integrating over space and using the Gauss Theorem, 

~J Gdv + J F-ds = 0 (2.2) 

V 3 

where v is any closed volume element and s is its surface. In the finite volume discretization, 
the computational domain of the problem is covered by a body conforming grid, typically 
of quadrilateral bricks. For the current scheme the grid need not be index separable. We 
only need a grid which provides reasonable resolution of the flow field. Finite element type 
of grids are acceptable as are also the well structured curvilinear grids used by many of 
the finite difference codes. Euler implicit discretization of the equation (2.2) is written as 

<w.« (- ’"‘j c " ) ^ (g»+‘) - D ixi (C"* 1 ) = 0 (2.3) 

where 

7j,k,t{G) = (F .S )j + i/2,k,i ~ S )> — 1/2 ,k,l ~ (F '.S ) }l k + l/2,t 

~{F .S )j t k-l/2,i + {F-s)j i k,t+i/2 - (F -S ) h k,i- 1/2 1 

Sj+i/2,k,i = (nds) ; a.x/ 2 ,k,i i etc - 

and n is the unit normal vector pointing in the increasing j direction. In the above, 
Dj,k,i{G) are the artificial dissipation terms similar to those used by Jameson, et al^ 5 ^. 
They are smaller than the trucation error terms of the scheme and are designed to suppress 
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the nonlinear instabilities which arise due to the central difference approximation of the 
convective flux terms, without affecting the accuracy of the solution. The implicit difference 
equations (2.3) for all the cells of the domain represent a set of coupled non-linear algebraic 
equations. The coupling of the equations for each cell is primarily with the equations of the 
neighboring nodes. Thus one can devise Gauss-Seidel type of iterations where equations 
at a group of neighboring nodes can be solved implicitly with the latest available iterates 
from the neighboring zones. The simplest scheme would be to solve the linearized form of 
the equations (2.3) one at a time with an inner iteration loop. 

Locally implicit approximation: The equation (2.3) is linearized and written in the form 

= ««”,*,< ( 2 - 4 ) 

where 

AG = G n+1 - G n (2.5) 

L(AG) ]iki t = CC • A G J)M 

+ CJM • &Gj-i t k,i + CJP • AGj+i t k,t 
+ D J M • AGj-2,k,i + DJ P • A Gj+2,k,t 
+ CKM * AGj t k-\,i + CKP • &Gj'k+i,i 
+ DKM • AGj t k-2,t + DKP ■ AG Jt k+2,e 
+ CLM • AG] t k,i-\ ■+■ CLP • AGj'k.t+i 
-I- DLM • AG h k,t-i -r DLP ■ AGj^^+i 

««?.*.« = -?i,kAG n ) + t>i,kAG n ) (2.7) 

The coefficients CJM.CJP, etc. are matrices which include the Jacobian matrices on the 
cell faces of the ( j,k,l ) cell and the dissipation coefficients. DJXf, DJP etc., are only 
diagonal matrices with artificial dissipation coefficients at the cell faces of {j,k, £), and CC 
is a matrix which can be approximated by a diagonal matrix which includes the time step 
term and the dissipation coefficients, both real and artificial. The linearized equations (2.4) 
are solved by an inner iteration for each time step by a modified symmetric Gauss-Seidel 
iteration. 

c ■ dG^l = R«” m - £(AG< >) J|M (2.8) 

A G (m) = A G (m “ l) + dG (m) , m = 1 8 (2.9) 

The superscript ( ) for AG inside the terms of the L operator of (2.8) is either (m) or 
(m — 1) whichever is the latest available value depending on the direction of the iteration 
sweep. The coefficient C is a modified form of the diagonal matrix CC at j, k , l designed 
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to enhance the stability and convergence characteristics of the scheme for large Courant 
numbers. Local time stepping enhances the convergence of the scheme for obtaining the 
steady state solutions. The eight inner iterations start from each of the eight corners of 
the computational domain to make the scheme symmetric. The initial A G^ can either be 
zero or can be the AG computed in the previous time step. The symmetric inner iteration 
is stable and removes any possible instabilities which may arise due to sweep direction 
being opposite to the flow direction. 

The equation (2.8) is now a set of scalar equations and we have explicit expressions 
for computing the corrections to the solution at each point. Thus the scheme is globally 
explicit while the inner iteration gives it a locally implicit character. Stability analysis of 
the model equations indicates that the scheme is unconditionally stable. In practice inner 
and outer relaxation parameters are introduced for equations (2.5) and (2.9) as follows: 

G n+X = G n + u 0Ut • AG (2.10) 

AG (m) = AG (m_1) + w tn • dG (m) (2.11) 

Under relaxation for w out in the range of 0.67 to 0.95 and over relaxation for u/ in from 1 
to 1.2 are indicated for good convergence of the scheme at Courant numbers of order 10 
and local time stepping. 


3. Multi-Grid Technique 

One of the most effective methods to speed up the convergence of iterative procedures 
for the solution of partial differential equations is the multi-grid technique. From a chosen 
fine grid, a sequence of coarser grids are formed by combining groups of adjacent cells 
in a systematic manner. For example, by omitting every other mesh line in an even cell 
numbered grid one can produce a coarse grid and repeating the process the next coarser 
level grid can be obtained. A set of iteration equations are solved on each of the grids 
where information from the fine grid is passed to the coarser grids and vice versa in each 
cycle. The highest frequency errors corresponding to each grid are dissipated rapidly in the 
iteration process at that grid level and thus errors of all frequencies are damped much more 
rapidly in the multi-grid mode than in a single grid iteration. For fluid flow problems, Ni^ 
and Jameson^ 7 ) have applied the multi-grid technique for inviscid transonic flow problems. 
A summary of the multi-grid technique is outlined below. 

At a particular grid level denoted by a mesh spacing parameter h, one or more time 
step calculations are earned out and the solution is updated. The solution G/i us injected 
to a coarser grid level indicated by the subscript 2h by the volume weighted average. 

G 2 h = v h G\) /V 2 k (3.1) 
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and the residual is injected by the summation over the fine grid cells, 


Rhk = ^ R-h 


A forcing function FF is defined on the coarse grid as the difference between the injected 
residual and the Res function (equation (2.7)) calculated on that grid using the injected 
solution. 

FF 2 h — Rhh ~ Res 2 h (3.3) 

The residual for advancing the solution by one or more time steps on the coarse grid is 
defined by 

R^h = Rz&2h + F F2h (3-4) 

The linearized equation corresponding to (2.4) on the coarse grid is 

L{AG) 2h = R 2tl (3.5) 

which is solved for one or more time steps by the locally implicit method outlined in the 
previous section. This process is repeated till the coarsest grid solution is completed. The 
coarse grid information is passed back to the finer grids as follows. The difference in the 
injected solution and the computed solution on the coarse grid is interpolated at the grid 
points of the next finer mesh and is added to the solution which was previously computed 
at that level. The solution is recomputed at that level for one or more time steps by 
the equations (3.4) and (3.5), using the forcing function FF which already exists at that 
level. This process is continued till the finest grid level is reached where the R and Res 
functions are the same. This completes one multi-grid cycle. In practice there are many 
variations to this method. One of the effective methods for fluid flows involves the solution 
advancement at each grid level for only one time step. Multi-grid scheme propagates the 
information from one boundary to another in each cycle and allows it to adjust to the 
boundary conditions quickly to reach a steady-state solution. 

4. Three-Dimensional Diffusion Equation 

As a test case for the performance of the algorithm, a time dependent diffusion equa- 
tion is solved with Dirichlet boundary conditions. 

d\L 

— = i/V 2 u, inD : 0 < x, y,z < 1 (4.1) 


u = 0 on dD 

At t = 0, u = Random numbers in (-1, 1) 
at all grid points 
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The solution u — * 0 as t —* oo. This problem has components at all frequencies supported 
by the grid and provides the most severe test for the process of convergence of the solution 
to steady state. The discretization chosen is the Euler implicit scheme, given by 


_L_ 

At 




v 

h? 


(* 


• +*, 1 


+ tf) 


(4.2) 


where A is forward difference and S 2 is central difference operator. The scheme has no 
limitation on grid sizes, which can be non-uniform in all directions. Constant grid size h 
is indicated here for the presentation of results for the test cases. The scheme (4.2) can be 
rewritten in the A-form as follows: 


{I-R(Sl + 61 + <J) } 1 = R'sl jtk (4.3) 

*«?•>,» = * («2 + + 

R h 2 

Locally implicit approximation : The implicit equations (4.3) are solved for Au’s by one 
or more Gauss-Seidel iterations. This test problem being relatively simple, one Gauss- 
Seidel iteration per time step without any further modification is adequate. It is stable 
for any R and is insensitive to the iteration sweep direction, though 8 symmetric sweeps 
(starting from each of the eight corners of the computational domain) per time step are 
more desirable for difficult problems, particularly with convection. The iteration for (4.3) 
can be written as 


(1 + 6R)Au\^ k = -r R (Au -r Auj j >Jfc _ 1 (4.4) 


^ U i+l,j,k + ^ u i,j+l,k ^ u ijU + l ) 


.( ) 


( ) 


where the superscript ( ) for Au on the right side is either (m) or (m - 1) depending 
on the iteration sweep direction. A simple iteration would be with m = l, Au^ 0 ^ = 0 
and the iteration sweep carried out in the directions of increasing indices. The multi-grid 
technique outlined in the last equation is applicable for this problem and is simpler since 
the equation (4.1) is scalar and linear. In each multi-grid cycle, the number of iterations 
carried out while going from fine to coarse grids are 1,2,2 2 , etc., and the same number 
in the reverse process. Figures 1-3 show the convergence of the solution u from random 
numbers to its asymptotic state u = 0. DELU = Max u — u aa y| = Mai jui. Work unit 
is the work required to compute the solution on the finest grid for one time step. For a 
20 3 mesh, convergence to more than 4 orders of magnitude is achieved in 60 work units 
for a single grid method where as it takes only 18 work units in a multi-grid mode with 
3 grids. With 40 3 mesh it takes 500 iterations for the same level of convergence with a 
single grid while it takes only 21 work units in the multi-grid mode with 4 grids. Even 
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with 40 x 40 x 80 mesh it takes only 21 work units with a 4 level multi-grid scheme while 
it takes more than 600 work units by a single grid scheme. Fig. 4 shows the convergence 
of Max \Residue\ against work units for a 40 3 mesh, which again shows a convergence of 
4 orders of magnitude in 21 work units. 


The locally implicit scheme can also be applied for a group of nodes at a time instead 
of one node at a time^ 5 ). For example, the equations (4.3) at a group of eight adjacent 
nodes can be solved with a locally implicit approximation. Fig. 5 shows the convergence 
of such a zonal scheme in both single grid and multi-grid modes. While the convergence 
is slightly better in terms of number of work units than the single point scheme, actual 
computational work is many times larger, since the locally implicit equations for each zone 
have to be solved simultaneously which is time consuming. Thus the single point scheme 
seems to be the computationally efficient scheme. 


5. Conclusions 

The locally implicit scheme is a computationally efficient scheme which converges 
rapidly in multi-grid modes for elliptic problems. It has the promise of providing a rapidly 
converging algorithm for steady-state viscous flow problems. 
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